//============================================================================
// Name        : sos1.cpp
// Author      : 
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C, Ansi-style
//============================================================================
#include <ilcplex/cplex.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define EPS 1e-6

int numcut1=0;

struct sosinfo {
	int 	numsos;		// how many SOS1 sets
	int		*sosbeg;	// length = numsos+1;
	int		*sosind;	// column index
	int		*sospos;	// sospos[column index] = position in sosind;
};

typedef struct sosinfo SOSINFO, *SOSINFOptr;

static int
	getImpSOS (CPXENVptr env, CPXLPptr lp, SOSINFOptr sosinfo);

static int CPXPUBLIC
	soscutscallback (CPXCENVptr env,
					 void 		*cbdata,
					 int 		wherefrom,
					 void 		*cbhandle,
					 int 		*useraction_p);

static int
	populateModel (CPXENVptr env, CPXLPptr lp);

int check=0;

int main (int  argc, char *argv[])
{
	int 	status = 0;

	CPXENVptr	env = NULL;
	CPXLPptr	lp  = NULL;

	SOSINFO sosinfo;

	double starttime, endtime, objval;

	env = CPXopenCPLEX (&status);
	status = CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_ON);

	lp = CPXcreateprob (env, &status, "myProb");
	populateModel(env, lp);
	status = CPXwriteprob (env, lp, "MIP.lp", NULL);

	status = CPXsetintparam (env, CPX_PARAM_PRELINEAR, 0);
	status = CPXsetintparam (env, CPX_PARAM_MIPCBREDLP, CPX_OFF);

//	status = CPXsetintparam (env, CPX_PARAM_HEURFREQ, -1);
//	status = CPXsetdblparam (env, CPX_PARAM_EPGAP, 0);
//	status = CPXsetintparam (env, CPX_PARAM_PREIND, 0);

	if ( status )  goto TERMINATE;

	status = CPXsetintparam (env, CPX_PARAM_THREADS, 1);
	status = CPXsetintparam (env, CPX_PARAM_MIPSEARCH, 1);

	status = getImpSOS (env, lp, &sosinfo);

	status = CPXsetcutcallbackfunc (env, soscutscallback, &sosinfo);

	status = CPXgettime(env, &starttime);

	if ( status )  goto TERMINATE;
	status = CPXmipopt (env, lp);

	status = CPXgetobjval(env, lp, &objval);
	printf("Objval: %f\n", objval);

	status = CPXgettime(env, &endtime);
	printf("Spend time: %f\n", endtime - starttime);

TERMINATE:
	if ( lp != NULL ) {
	   status = CPXfreeprob (env, &lp);
	   if ( status ) {
		  fprintf (stderr, "CPXfreeprob failed, error code %d.\n",
				   status);
	   }
	}

	if ( env != NULL ) {
	   status = CPXcloseCPLEX (&env);
	   if ( status ) {
		  char errmsg[1024];
		  fprintf (stderr, "Could not close CPLEX environment.\n");
		  CPXgeterrorstring (env, status, errmsg);
		  fprintf (stderr, "%s", errmsg);
	   }
	}

	return (status);
}

static int
	getImpSOS(CPXENVptr env, CPXLPptr lp, SOSINFOptr sosinfo)
{
	int status = 0;
	int numrows, numcols, rnzcnt, rmatspace, surplus;

	int 	*sosstatus 	= NULL;

	int 	*rmatind 	= NULL,
			*rmatbeg 	= NULL;
	double 	*rmatval 	= NULL,
			*rhs 		= NULL;
	char 	*vType   	= NULL;
	char	*sense	 	= NULL;

	int ifcardcon;

	int i, j, numsos;

	numcols = CPXgetnumcols(env, lp);
	numrows = CPXgetnumrows(env, lp);

	status  = CPXgetrows(env, lp, &rnzcnt, NULL, NULL, NULL, 0, &surplus, 0, numrows-1);
	rmatspace = -surplus;
	rmatbeg = (int*) malloc(numrows*sizeof(int));
	rmatind = (int*) malloc(rmatspace*sizeof(int));
	rmatval = (double*) malloc(rmatspace*sizeof(double));
	status  = CPXgetrows(env, lp, &rnzcnt, rmatbeg, rmatind, rmatval, rmatspace, &surplus, 0, numrows-1);

	rhs = (double*) malloc(numrows*sizeof(double));
	status = CPXgetrhs(env, lp, rhs, 0, numrows-1);

	sense = (char*) malloc(numrows*sizeof(char));
	status = CPXgetsense(env, lp, sense, 0, numrows-1);

	vType  = (char*) malloc (numcols*sizeof(char));
	status = CPXgetctype(env, lp, vType, 0, numcols-1);

	sosstatus = (int*) malloc(numcols*sizeof(int));
	for (i=0; i<numcols; i++) sosstatus[i]=-1;

	numsos=0;
	for (i=0; i<numrows; i++) {
		if ((abs(rhs[i]-1)<=EPS && sense[i]=='L') || (abs(rhs[i]+1)<=EPS && sense[i]=='G')) {
			ifcardcon=1;
			for (j=rmatbeg[i]; j < (i+1 < numrows? rmatbeg[i+1]: rnzcnt); j++) {
				if (vType[rmatind[j]]!='B') {ifcardcon=0; break;}
				if (sense[i]=='L' && abs(rmatval[j]-1)>=EPS) {ifcardcon=0; break;}
				if (sense[i]=='G' && abs(rmatval[j]+1)>=EPS) {ifcardcon=0; break;}
			}
			if (ifcardcon) {
				numsos++;
			}
		}
		// link the variable bound of x <= y, such that sosstatus[colind[y]] = colind[x];
		if ((i+1 < numrows? rmatbeg[i+1]: rnzcnt) - rmatbeg[i]==2 && abs(rhs[i])<=EPS) {
			j = rmatbeg[i];
			if (sense[i] == 'L' &&
				vType[rmatind[j]]=='B' && abs(rmatval[j]+1)<=EPS && vType[rmatind[j+1]]=='C' && abs(rmatval[j+1]-1)<=EPS)
			{
				sosstatus[rmatind[j]] = rmatind[j+1];
			}
			if (sense[i] == 'L' &&
				vType[rmatind[j]]=='C' && abs(rmatval[j]-1)<=EPS && vType[rmatind[j+1]]=='B' && abs(rmatval[j+1]+1)<=EPS)
			{
				sosstatus[rmatind[j+1]] = rmatind[j];
			}
			if (sense[i] == 'G' &&
				vType[rmatind[j]]=='B' && abs(rmatval[j]-1)<=EPS && vType[rmatind[j+1]]=='C' && abs(rmatval[j+1]+1)<=EPS)
			{
				sosstatus[rmatind[j]] = rmatind[j+1];
			}
			if (sense[i] == 'G' &&
				vType[rmatind[j]]=='C' && abs(rmatval[j]+1)<=EPS && vType[rmatind[j+1]]=='B' && abs(rmatval[j+1]-1)<=EPS)
			{
				sosstatus[rmatind[j+1]] = rmatind[j];
			}
		}
	}

	sosinfo->sosbeg = (int*) malloc((numsos+1)*sizeof(int));
	sosinfo->sosind = (int*) malloc(numcols*sizeof(int));
	sosinfo->sospos = (int*) malloc(numcols*sizeof(int));
	for (j=0; j<numcols; j++) {
		sosinfo->sospos[j] = -1;
	}

	numsos = 0;
	sosinfo->sosbeg[0]=0;
	for (i=0; i<numrows; i++) {
		if ((abs(rhs[i]-1)<=EPS && sense[i]=='L') || (abs(rhs[i]+1)<=EPS && sense[i]=='G')) {
			ifcardcon=1;
			for (j=rmatbeg[i]; j < (i+1 < numrows? rmatbeg[i+1]: rnzcnt); j++) {
				if (vType[rmatind[j]]!='B') {ifcardcon=0; break;}
				if (sense[i]=='L' && abs(rmatval[j]-1)>=EPS) {ifcardcon=0; break;}
				if (sense[i]=='G' && abs(rmatval[j]+1)>=EPS) {ifcardcon=0; break;}
			}
			if (ifcardcon) {
				sosinfo->sosbeg[numsos+1] = sosinfo->sosbeg[numsos];
				for (j=rmatbeg[i]; j < (i+1 < numrows? rmatbeg[i+1]: rnzcnt); j++) {
					sosinfo->sosind[sosinfo->sosbeg[numsos+1]] = sosstatus[rmatind[j]];
					sosinfo->sospos[sosstatus[rmatind[j]]] = sosinfo->sosbeg[numsos+1];
					sosinfo->sosbeg[numsos+1]++;
				}
				numsos++;
			}
		}
	}

	sosinfo->numsos=numsos;

//	for (i=0; i<numsos; i++) {
//		for (j=sosinfo->sosbeg[i]; j<sosinfo->sosbeg[i+1]; j++) {
//			printf("col = %d ", sosinfo->sosind[j]);
//		}
//		printf("\n");
//	}
TERMINATE:
	free(sosstatus);

	free(rmatbeg);
	free(rmatind);
	free(rmatval);
	free(rhs);
	free(sense);
	free(vType);

	return status;
}

static int CPXPUBLIC
soscutscallback (CPXCENVptr env,
               	 void       *cbdata,
               	 int        wherefrom,
               	 void       *cbhandle,
               	 int        *useraction_p)
{
	int status  = 0;

	SOSINFOptr sosinfo = (SOSINFOptr) cbhandle;

	int numrows, numcols, rnzcnt, rmatspace, surplus;

	char	**colname = NULL,
			**rowname = NULL;
	char	*colnamestore = NULL,
			*rownamestore = NULL;
	int		colnamespace, rownamespace;

	int		numsos 	= sosinfo->numsos;
	int		*sosbeg	= sosinfo->sosbeg,
			*sosind = NULL,
			*sospos = sosinfo->sospos;

	int 	*rmatind = NULL,
			*rmatbeg = NULL,
			*cutind	 = NULL;
	double 	*x		 = NULL,
			*rmatval = NULL,
			*cutval  = NULL,
			*soswt 	 = NULL,
			*rhs 	 = NULL;
	char	*sense	 = NULL;

	int i,j,k,l,col,tempind,sosrow,cutnz;

	double tempval,cutrhs,cutvio,tempsum,a1;

	int sossize = sosinfo->sosbeg[numsos];

	*useraction_p = CPX_CALLBACK_DEFAULT;

	CPXCLPptr lp;
	status = CPXgetcallbacklp (env, cbdata, wherefrom, &lp);

	numcols = CPXgetnumcols(env, lp);
	numrows = CPXgetnumrows(env, lp);

	status = CPXgetcolname (env, lp, NULL, NULL, 0, &surplus, 0, numcols-1);
	colnamespace = - surplus;
	colname      = (char **) malloc (sizeof(char *)*numcols);
	colnamestore = (char *)  malloc (colnamespace);
	status = CPXgetcolname (env, lp, colname, colnamestore, colnamespace, &surplus, 0, numcols-1);

	status  = CPXgetrows(env, lp, &rnzcnt, NULL, NULL, NULL, 0, &surplus, 0, numrows-1);
	rmatspace = -surplus;
	rmatbeg = (int*) malloc(numrows*sizeof(int));
	rmatind = (int*) malloc(rmatspace*sizeof(int));
	rmatval = (double*) malloc(rmatspace*sizeof(double));
	status  = CPXgetrows(env, lp, &rnzcnt, rmatbeg, rmatind, rmatval, rmatspace, &surplus, 0, numrows-1);

	rhs = (double*) malloc(numrows*sizeof(double));
	status = CPXgetrhs(env, lp, rhs, 0, numrows-1);

	sense = (char*) malloc(numrows*sizeof(char));
	status = CPXgetsense(env, lp, sense, 0, numrows-1);

	x  = (double*) malloc (numcols*sizeof(double));
	status = CPXgetcallbacknodex  (env, cbdata, wherefrom, x, 0, numcols-1);

	sosind = (int*) malloc(sossize*sizeof(int));
	soswt  = (double*) malloc(sossize*sizeof(double));

	cutind = (int*) malloc(numcols*sizeof(int));
	cutval = (double*) malloc(numcols*sizeof(double));

	for (i=0; i<numrows; i++) {
		// TODO need more theories to add cuts on more complicated rows
		sosrow = 1;
		for (j=rmatbeg[i]; j < (i+1 < numrows? rmatbeg[i+1]: rnzcnt); j++) {
			if ( sospos[rmatind[j]] == -1 ||
				 (sense[i] == 'L' || sense[i] == 'E') && rmatval[j] < 0 ||
				 (sense[i] == 'G' || sense[i] == 'E') && rmatval[j] > 0 ) {
				sosrow = 0;
				break;
			}
		}
		if (!sosrow) continue;

		for (j=0; j<sossize; j++) soswt[j] = 0.0;
		for (j=rmatbeg[i]; j < (i+1 < numrows? rmatbeg[i+1]: rnzcnt); j++) {
			col = rmatind[j];
			if (sospos[col] != -1) soswt[sospos[col]] = rmatval[j];
		}

		for (j=0; j<sossize; j++) sosind[j] = sosinfo->sosind[j];

		for (j=0; j<numsos; j++) {
			for (k=sosbeg[j]; k<sosbeg[j+1]; k++) {
				for (l=k+1; l<sosbeg[j+1]; l++) {
					if (soswt[k] < soswt[l]) {
						tempind = sosind[k];
						sosind[k] = sosind[l];
						sosind[l] = tempind;
						tempval = soswt[k];
						soswt[k] = soswt[l];
						soswt[l] = tempval;
					}
				}
			}
		}

//		if (check==0) {
//		for (j=0; j<numsos; j++) {
//			for (k=sosbeg[j]; k<sosbeg[j+1]; k++) {
//				printf("col = %d val = %f ", sosind[k], soswt[k]);
//			}
//			printf("\n");
//		}
//		check=1;
//		}

		//add cuts from theory 1.
		tempsum = 0.0;
		for (j=0; j<numsos; j++) {
			for (k=sosbeg[j]; k<sosbeg[j+1]; k++) {
				if (abs(x[sosind[k]]) >= EPS && soswt[k] > 0) {
					tempsum += soswt[sosbeg[j]];
					break;
				}
			}
		}

		if (tempsum <= rhs[i]) continue;

		cutrhs = rhs[i];
		cutvio = -cutrhs;
		cutnz = 0;
		for (j=0; j<numsos; j++) {
			a1 = 0;
			for (k=sosbeg[j]; k<sosbeg[j+1]; k++) {
				if (abs(x[sosind[k]]) >= EPS && soswt[k] > 0) {
					a1=soswt[k];
					break;
				}
			}
			for (l=k; l<sosbeg[j+1];l++) {
				if (soswt[l] > 0) {
					cutind[cutnz] = sosind[l];
					cutval[cutnz] = (soswt[l] > rhs[i] - tempsum + a1)? soswt[l]: rhs[i] - tempsum + a1;
					cutvio += soswt[cutnz]*x[cutind[cutnz]];
					cutnz++;
				}
			}
		}

		if (cutvio > 0.01 && numcut1 < 1000) {
			status = CPXcutcallbackadd (env, cbdata, wherefrom, cutnz, cutrhs, 'L', cutind, cutval, 0);
			if ( status ) {
				fprintf (stderr, "Failed to add cut.\n");
				goto TERMINATE;
			}
//			numcut1++;
//			if (numcut1==2) {
//				printf("i=%d\n",i);
//				for (j=0; j<cutnz; j++) {
//					printf("%f %s +", cutval[j],colname[cutind[j]]);
//				}
//				printf("\n");
//			}
		}
	}

	*useraction_p = CPX_CALLBACK_SET;

TERMINATE:
	free(rmatbeg);
	free(rmatind);
	free(rmatval);
	free(rhs);
	free(sense);

	free(cutind);
	free(cutval);

	return (status);
}




static int
	populateModel (CPXENVptr env, CPXLPptr lp)
{
	int status;
	int numSOS  = 10,
		numrows = 100,
		soslb = 5,
		sosub = 10,
		density = 20;

	int *numSize, *rmatind;
	int **xcol, **ycol, **sosobj, **a;
	int i,j,k,numcols=0,nzcnt,zero=0;

	double l=0, u=1, *rmatval, rhs, max_a, colobj;
	char *colname, *rowname, sense='L', vch='C';
	srand(1);

	numSize = (int*) malloc(numSOS*sizeof(int));
	// each SOS set has 10 ~ 20 members
	for (i=0; i<numSOS; i++) {
		numSize[i] = rand() % (sosub-soslb+1) + soslb;
	}

	xcol = (int**) malloc(numSOS*sizeof(int*));
	ycol = (int**) malloc(numSOS*sizeof(int*));

	CPXchgobjsen (env, lp, CPX_MAX);

	sosobj = (int**) malloc(numSOS*sizeof(int*));
	for (i=0; i<numSOS; i++) {
		sosobj[i] = (int*) malloc(numSize[i]*sizeof(int));
		for (j=numSize[i]-1; j>=0; j--) {
			if (j==numSize[i]-1) sosobj[i][j] = rand()%4+1;
			else sosobj[i][j] = sosobj[i][j+1] + rand()%4+1;
		}
	}
	
	colname = (char*) malloc(30*sizeof(char));
	for (i=0; i<numSOS; i++) {
		xcol[i] = (int*) malloc(numSize[i]*sizeof(int));
		ycol[i] = (int*) malloc(numSize[i]*sizeof(int));
		for (j=0; j<numSize[i]; j++) {
			sprintf(colname,"x_%d_%d",i,j);
			colobj = sosobj[i][j];
			vch='C';
			status = CPXnewcols(env, lp, 1, &colobj, &l, &u, &vch, &colname);
			xcol[i][j]=numcols;
			numcols++;
			sprintf(colname,"y_%d_%d",i,j);
			colobj = 0;
			vch='B';
			status = CPXnewcols(env, lp, 1, &colobj, &l, &u, &vch, &colname);
			ycol[i][j]=numcols;
			numcols++;
		}
	}
	free(colname);

	a = (int**) malloc(numSOS*sizeof(int*));
	for (i=0; i<numSOS; i++) {
		a[i] = (int*) malloc(numSize[i]*sizeof(int));
	}

	rowname = (char*) 	malloc(30*sizeof(char));
	rmatind = (int*)  	malloc(numcols*sizeof(int));
	rmatval = (double*) malloc(numcols*sizeof(double));

	for (i=0; i<numrows; i++) {
		sprintf(rowname,"c_%d",i);
		nzcnt = 0;
		max_a = 0;
		rhs = 0.0;

		for (j=0; j<numSOS; j++) {
			for (k=numSize[j]-1; k>=0; k--) {
				if (k==numSize[j]-1) a[j][k] = rand()%3+1;
				else a[j][k] = ceil(sosobj[j][k]*a[j][k+1]/sosobj[j][k+1]) + 1 + rand()%5;
			}
			if (rand()%100>90) continue;
			for (k=0; k<numSize[j]; k++) {
				if (rand()%100>density) continue;
				rmatind[nzcnt] = xcol[j][k];
				rmatval[nzcnt] = a[j][k];
				if (max_a < rmatval[nzcnt]) max_a = rmatval[nzcnt];
				nzcnt++;
			}
			rhs += a[j][numSize[j]-1];
		}
		rhs = floor(rhs * 1.3);
		if (rhs < max_a + 1) rhs = max_a + 1;
		status = CPXaddrows(env, lp, 0, 1, nzcnt, &rhs, &sense, &zero, rmatind, rmatval, NULL, &rowname);
	}

	for (i=0; i<numSOS; i++) {
		for (j=0; j<numSize[i]; j++) {
			sprintf(rowname,"bound_%d_%d",i,j);
			rmatind[0] = xcol[i][j]; rmatind[1] = ycol[i][j];
			rmatval[0] = 1; rmatval[1] = -1;
			rhs = 0.0;
			status = CPXaddrows(env, lp, 0, 1, 2, &rhs, &sense, &zero, rmatind, rmatval, NULL, &rowname);
		}
	}

	for (i=0; i<numSOS; i++) {
		sprintf(rowname,"card_%d",i);
		for (j=0; j<numSize[i]; j++) {
			rmatind[j] = ycol[i][j];
			rmatval[j] = 1;
		}
		rhs=1.0;
		status = CPXaddrows(env, lp, 0, 1, numSize[i], &rhs, &sense, &zero, rmatind, rmatval, NULL, &rowname);
	}

	free(rowname);
	free(rmatind);
	free(rmatval);
	free(a);
	free(xcol);
	free(ycol);
	free(sosobj);

	return status;
}
